#still working on this file

getAllSNPs <- function(embl.file,grouseMapperDataFrames,tsvFiles,tsvLineDataFrames){
	#embl.file contains basic annotation information, including what the contigs are
	#grouseMapperDataFrames is a list of dataframes produced by readNewblerMapper
	#tsv files is a character vector containing names of tsv files
	#tsvLineDataFrames is a list of dataframes produced by readTsvFile
	
	embl.lines <- readLines(embl.file)
	embl.CDS.lines <- grep('^FT[ ]*CDS',embl.lines)
	embl.cDNA <- embl.lines[grep('^ID',embl.lines)]
	embl.cDNA <- sub('^ID[ ]*','',embl.cDNA)
	embl.cDNA <- sub('_matched$','',embl.cDNA)
	
	embl.product.lines <- grep('/product',embl.lines)
	embl.CDS <- rep('xx',length(embl.product.lines))
	for(cdsi in 1:length(embl.product.lines)){
		tmp <- embl.lines[embl.CDS.lines[cdsi]:(embl.product.lines[cdsi]-1)]
		tmp <- sub('^FT[ ]*','',tmp)
		tmp <- sub('^CDS[ ]*','',tmp)
		embl.CDS[cdsi] <- paste(tmp,collapse="")
		}
	#embl.CDS <- gsub('FT[ ]*','',embl.CDS)
	#embl.CDS <- sub('CDS[ ]*','',embl.CDS)
	embl.CDS.dir <- "forward"
	if(substr(embl.CDS[1],1,4)=="comp") embl.CDS.dir <- "reverse"
	embl.CDS <-sub('[a-z]*\\(','',embl.CDS)
	embl.CDS <- sub('\\)$','',embl.CDS)

	#create a complex list
	embl.CDS.list <- lapply(strsplit(embl.CDS,','),strsplit,'\\.\\.')
	embl.CDS.start <- lapply(embl.CDS.list,function(X){sapply(X,function(X){as.numeric(X[1])})})
	embl.CDS.stop <- lapply(embl.CDS.list,function(X){sapply(X,function(X){as.numeric(X[2])})})

	embl.length <- as.numeric(sub('^FT[ ]*source[ ]*1\\.\\.','',embl.lines[grep('^FT[ ]*source',embl.lines)]))

	
	contig.names <- strsplit(paste(sub('^CC[ ]*','',embl.lines[grep('^CC',embl.lines)]),collapse=" "),'; ')[[1]]

	#where are the contigs
	contig.lines <- grep('/note=\"contig',embl.lines)
	contig.direction <- 	factor(rep("forward",length(contig.lines)),levels=c("forward","reverse"))
	contig.dir.txt <- embl.lines[contig.lines+1]
	contig.direction[grep('reverse',contig.dir.txt)] <- "reverse"
	
	contig.pos <- sub('^FT[ ]*misc_feature[ ]*','',embl.lines[contig.lines-1])
	contig.mat <- sapply(contig.pos,function(X){as.numeric(strsplit(X,'\\.\\.')[[1]])})
	#1st row has start points, 2nd row has end points

	#generate a big list of coverage & snps
	contig.tsv <- vector('list',length=length(grouseMapperDataFrames))
	contig.snp <- vector('list',length=length(grouseMapperDataFrames))
	names(contig.tsv) <- paste('pop',1:length(contig.tsv),sep="")
	names(contig.snp) <- paste('pop',1:length(contig.snp),sep="")

	#first start by collecting snp and coverage info

	for(pop in 1:length(contig.tsv)){

		for(cnti in 1:length(contig.names)){
			contig.tsv[[pop]][[cnti]] <- getContigFromTsv(contig.names[cnti], tsvLineDataFrames[[pop]],tsvFiles[pop])
			
			
			names(contig.tsv[[pop]])[cnti] <- names(contig.tsv)[cnti]
			
			#trap errors when empty data.frame returned
			tmp.try <- grouseMapperDataFrames[[pop]][grouseMapperDataFrames[[pop]]$ref.accno==contig.names[cnti],]
			if(nrow(tmp.try)==0){
				tmp.try$tsv.coverage <- numeric(0)
				tmp.try$rel.pos <- numeric(0)
				tmp.try$freq <- numeric(0)
				contig.snp[[pop]][[cnti]] <- tmp.try
				next # JUMP TO NEXT cnti
				}
			
			contig.snp[[pop]][[cnti]] <- tmp.try
			contig.snp[[pop]][[cnti]]$tsv.coverage <- 0
			
			#this line picks coverage out of the contig.tsv df and adds it to contig.snp
			#check carefully!
			#contig.snp[[pop]][[cnti]]$tsv.coverage[ contig.snp[[pop]][[cnti]]$start.pos %in%  contig.tsv[[pop]][[cnti]]$position] <- contig.tsv[[pop]][[cnti]]$depth[ match(contig.snp[[pop]][[cnti]]$start.pos, contig.tsv[[pop]][[cnti]]$position) ]
	
			if(contig.direction[cnti]=="forward"){
				contig.snp[[pop]][[cnti]]$rel.pos <- contig.mat[1, cnti] -1 + contig.snp[[pop]][[cnti]]$start.pos
				contig.tsv[[pop]][[cnti]]$rel.pos <- contig.mat[1, cnti] -1 + contig.tsv[[pop]][[cnti]]$position
			}else{ #reverse
				contig.snp[[pop]][[cnti]]$rel.pos <- contig.mat[2, cnti] - contig.snp[[pop]][[cnti]]$start.pos +1
				contig.tsv[[pop]][[cnti]]$rel.pos <- contig.mat[2, cnti] - contig.tsv[[pop]][[cnti]]$position +1
				}#end of if for/rev condition
			contig.snp[[pop]][[cnti]]$freq <- as.numeric(sub('%$','',contig.snp[[pop]][[cnti]]$freq))/100
			}#end of cnti loop

		names(contig.tsv[[pop]]) <- contig.names
		names(contig.snp[[pop]]) <- contig.names
		}#end of pop loop




	#collect all the SNPs together
	contig.snp2 <- vector('list',length=length(grouseMapperDataFrames))
	for(i in 1:length(contig.snp)){
		for(j in 1:length(contig.snp[[i]])){
			if(j==1){
				tmp <- contig.snp[[i]][[j]]
				}else{
				tmp <- rbind(tmp,contig.snp[[i]][[j]])
				}
			#
			
			}#end of j loop
		names(tmp)[names(tmp) %in% c("start.pos", "rel.pos")] <- c("old.pos","start.pos")
		
		if(nrow(tmp)==0){ #trap errors where nothing sequenced for a population
			tmp$cDNA <- character(0)
			tmp$id <- character(0)
			contig.snp2[[i]] <- tmp
			}else{
			tmp$cDNA <- embl.cDNA #crashes if tmp has zero rows
			tmp$id <- paste(tmp$cDNA,tmp$start.pos,tmp$ref.accno,tmp$old.pos,tmp$ref.nucl,tmp$var.nucl,sep="_")
			contig.snp2[[i]] <- tmp[order(tmp$start.pos),]
			}
		
		if(i == 1){ #collect together snp ids
			tmp.id <- as.character(contig.snp2[[i]]$id)
			}else{
			tmp.id <- c(tmp.id,as.character(contig.snp2[[i]]$id))
			}
		}
	
	snp.df1 <- data.frame(id = unique(tmp.id),stringsAsFactors=FALSE)
	snp.df1$cDNA <- embl.cDNA
	snp.df1$start.pos <- sapply(X=snp.df1$id,FUN=function(X){
		tmp <- strsplit(X,"_")[[1]]
		as.numeric(tmp[length(tmp)-4])
		})

	snp.df1$contig <- sapply(X=snp.df1$id,FUN=function(X){
		tmp <- strsplit(X,"_")[[1]]
		tmp[grep('contig',tmp)]
		})
	snp.df1$old.pos <- sapply(X=snp.df1$id,FUN=function(X){
		tmp <- strsplit(X,"_")[[1]]
		as.numeric(tmp[length(tmp)-2])
		})
	
	snp.df1 <- snp.df1[order(snp.df1$start.pos),]
	
	for(pop in 1:length(grouseMapperDataFrames)){
		tmp.nms <- c(names(snp.df1),paste("pop",pop,c("freq","mapper.depth","tsv.depth","quality.score","signal","std.deviation"),sep="."))
		snp.df1[,ncol(snp.df1)+1:6] <- 0
		names(snp.df1) <- tmp.nms
		}
	
	#associated SNPs with tsv.coverage data for each snp
	contigs.with.snps <- unique(snp.df1$contig)
	for(pop in 1:length(contig.snp)){
		for(cnti in contigs.with.snps){
			#get snp freq and mapper depth from contig.snp
			tmp.snp.df <- contig.snp[[pop]][[cnti]]
			tmp.tsv.df <- contig.tsv[[pop]][[cnti]]
			
			snp.df1[snp.df1$contig==cnti, match(paste("pop",pop,c("freq","mapper.depth"),sep="."),names(snp.df1))] <- 
				as.data.frame(t(sapply(X=snp.df1$old.pos[snp.df1$contig==cnti],FUN= function(X,tmp.snp.df){findSnpFreq(X,tmp.snp.df)},tmp.snp.df)))
			
			snp.df1[snp.df1$contig==cnti, match(paste("pop",pop,c("tsv.depth","quality.score","signal","std.deviation"),sep="."),names(snp.df1))] <- 
				as.data.frame(t(sapply(X=snp.df1$old.pos[snp.df1$contig==cnti],FUN= function(X,tmp.tsv.df){findTsvCov(X,tmp.tsv.df)},tmp.tsv.df)))

			}
		
		
		}

	snp.df1 
	}

#get freq and mapping coverage from contig.snp
findSnpFreq <- function(old.pos,contig.snp.df){
	if(class(contig.snp.df)!="data.frame") return(c(0,0))
	if(nrow(contig.snp.df)==0) return(c(0,0))
	if(!any(contig.snp.df$start.pos %in% old.pos)) return(c(0,0))
	as.numeric(contig.snp.df[match(old.pos,contig.snp.df$start.pos),c("freq","depth")])
	}

findTsvCov <- function(old.pos,contig.tsv.df){
	if(class(contig.tsv.df)!="data.frame") return(c(0,0,0,0))
	if(nrow(contig.tsv.df)==0) return(c(0,0,0,0))
	if(!any(contig.tsv.df$position %in% old.pos)) return(c(0,0,0,0))
	as.numeric(contig.tsv.df[match(old.pos,contig.tsv.df$position),c("depth","quality.score","signal","std.deviation")])
	}



#create a little function to find info for each snp
#findTsvCoverage <- function(contig.tsv,pop,contig,old.pos){
	#contig.tsv is a complex data object created internally by getAllSNPs
#	tmp.df <- contig.tsv[[pop]][[contig]]
#	if(class(tmp.df)!="data.frame") return(0)
#	tst.out <- tmp.df$depth[tmp.df$position==old.pos]
#	if(length(tst.out)==0){
#		warning("zero length returned for pop ",pop,", ",contig,", position ",old.pos,"\n")
#		return(0)
#		}
#	if(length(tst.out)>1){
#		warning("length >1 returned for pop ",pop,", ",contig,", position ",old.pos,"\n")
#		return(0)
#		}
#	as.numeric(tst.out)
#	}

#debug(getAllSNPs)
tst <- getAllSNPs(embl.file="/Paterson/Datafiles/grouse/embl_exon2/cDNA_469-1_exon.embl",grouseMapperDataFrames=list(grouse1= grouse1All,grouse2 = grouse2All,grouse3= grouse3All),tsvFiles=c('mapping/grouse1.tsv','mapping/grouse2.tsv','mapping/grouse3.tsv'),tsvLineDataFrames=list(grouse1=grouse1_lines,grouse2=grouse2_lines,grouse3=grouse3_lines))
tst2 <- codingSNPref(tst,seq.name="cDNA_469-1") #work out dn and ds













captureSNPs2 <- function(contig.map,contig.dataframes,cDNA.name,AceCOs,ace.files,envir=.GlobalEnv){
	#a wrapper to take SNP data from multiple populations
	#contig map is a list derived from join contigs
	#contig.dataframes is a character vector containing names of dataframes
	# produced by readNewblerMapper
	#cDNA.name is the name of the cDNA-derived supercontig
	#AceCOs is a character vector containing names of dataframes for CO line
	#numbers in ace files produced by readAceCO
	#ace.files is a character vector containing names of ace files to be read
	#by readContigAce
	#envir is the environment in which to look for contig.dataframes & AceCOs
	
	SNPlist <- list()
	AllCoverage <- list()
	
	for(i in 1:length(contig.dataframes)){
		contig.df <- get(contig.dataframes[i],envir=envir,inherits=FALSE)
		SNPlist[[i]] <- mapSNPsToContigs(contig.df,contig.map,cDNA.name)
		if(!identical(SNPlist[[i]],0)){
			SNPlist[[i]]$SNPid <- paste(cDNA.name,SNPlist[[i]]$start.pos,SNPlist[[i]]$ref.accno,SNPlist[[i]]$old.start,SNPlist[[i]]$ref.nucl,SNPlist[[i]]$var.nucl,sep="_")
			}else{
			SNPlist[[i]]$SNPid <- NULL
			}
		
		AceCO <- get(AceCOs[i],envir=envir,inherits=FALSE)
		ace.file<-ace.files[i]
		linesAceList <- lapply(X=contig.map[[cDNA.name]]$sc.contig,
			FUN=function(X,ace.file,AceCO){
				readContigAce(ace.file,AceCO,X)
				},ace.file,AceCO)
		AceObjectList <- lapply(X= linesAceList,FUN=makeAceObject)
		CoverageList <- lapply(X= AceObjectList,FUN= mapCoverage)
		rm(linesAceList, AceObjectList)
		names(CoverageList) <- contig.map[[cDNA.name]]$sc.contig
		#should now have an estimate of the coverage at each base
		AllCoverage[[i]] <- CoverageList
		
		}
	names(SNPlist) <- contig.dataframes
	names(AllCoverage) <- contig.dataframes
	AllSNPids <- sort(unique(unlist(sapply(X=SNPlist,FUN=function(X){
		X$SNPid
		}))))
	list(SNPlist = SNPlist ,AllCoverage=AllCoverage,AllSNPids=AllSNPids)
	}
